ENMTools

This package implements various tests, visualizations, and metrics for use with environmental niche models (ENMs) and species distribution models (SDMs).


Installation

At present, ENMTools is downloadable from https://github.com/danlwarren/ENMTools. There are multiple ways to download it. The easiest is to use devtools and install from GitHub.

Installing from GitHub using devtools

Run the following code from your R console:

install.packages("devtools")
library(devtools)
install_github("danlwarren/ENMTools")
library(ENMTools)

Install from zip file

A zipped version of the package is available at https://github.com/danlwarren/ENMTools/archive/master.zip. To install from the zip file, download a copy of it to your system. Once it’s finished downloading, type the following (where PATH is the path to the zip file):

install.packages("devtools")
library(devtools)
install_local("PATH")
library(ENMTools)

Interacting with ENMTools

Creating enmtools.species objects

First we’re going to load in some environmental data.

env.files <- list.files(path = "test/testdata/", pattern = "pc", full.names = TRUE)
env <- stack(env.files)
names(env) <- c("layer.1", "layer.2", "layer.3", "layer.4")
env <- setMinMax(env)

ENMTools is primarily designed to examine patterns of similarity and difference between ENMs for different species. In order to simplify interactions with the functions in ENMTools, you need to put your data for each of your species into an enmtools.species object. You can create and view an empty enmtools.species object just by typing:

ahli <- enmtools.species()
ahli
## 
## 
## Range raster not defined.
## 
## Presence points not defined.
## 
## Background points not defined.
## 
## Species name not defined.

You can add data to this object manually:

names(ahli)
## [1] "range"             "presence.points"   "background.points"
## [4] "models"            "species.name"
ahli$species.name <- "ahli"
ahli$presence.points <- read.csv("test/testdata/ahli.csv")[,3:4]
ahli$range <- background.raster.buffer(ahli$presence.points, 50000, mask = env)
ahli$background.points <- background.points.buffer(points = ahli$presence.points,
                                                   radius = 20000, n = 1000, mask = env[[1]])

ahli
## 
## 
## Range raster: 
## class       : RasterLayer 
## dimensions  : 418, 1535, 641630  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -86.90809, -74.11642, 19.80837, 23.2917  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer.1 
## values      : 1, 1  (min, max)
## 
## 
## 
## Presence points (first ten only): 
## 
##  Longitude   Latitude
## ----------  ---------
##   -80.0106    21.8744
##   -79.9086    21.8095
##   -79.8065    21.7631
##   -79.8251    21.8095
##   -79.8807    21.8374
##   -79.9550    21.8374
##   -80.3446    22.0136
##   -80.2983    21.9951
##   -80.1776    21.9023
##   -80.1591    21.9673
## 
## 
## Background points (first ten only): 
## 
##  Longitude   Latitude
## ----------  ---------
##  -80.43726   22.17087
##  -79.99559   22.02920
##  -79.89559   21.76254
##  -79.76226   21.92920
##  -80.22059   22.07087
##  -80.36226   22.01254
##  -80.01226   21.91254
##  -80.48726   22.11254
##  -80.07892   21.86254
##  -80.32059   21.96254
## 
## 
## Species name:  ahli

Or you can add bits of it when the object is created:

allogus <- enmtools.species(species.name = "allogus", 
                            presence.points = read.csv("test/testdata/allogus.csv")[,3:4])

allogus$range <- background.raster.buffer(allogus$presence.points, 50000, mask = env)
allogus$background.points <- background.points.buffer(points = allogus$presence.points,
                                                      radius = 20000, n = 1000, mask = env[[1]])

allogus
## 
## 
## Range raster: 
## class       : RasterLayer 
## dimensions  : 418, 1535, 641630  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -86.90809, -74.11642, 19.80837, 23.2917  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer.1 
## values      : 1, 1  (min, max)
## 
## 
## 
## Presence points (first ten only): 
## 
##  Longitude   Latitude
## ----------  ---------
##   -79.2527    22.2109
##   -78.7774    22.2241
##   -78.6189    22.2373
##   -78.1039    21.1809
##   -78.0247    21.1809
##   -77.9983    20.9301
##   -77.9719    21.7091
##   -77.9719    21.5507
##   -77.9323    21.6167
##   -77.9323    20.7320
## 
## 
## Background points (first ten only): 
## 
##  Longitude   Latitude
## ----------  ---------
##  -75.28726   20.02920
##  -75.86226   20.18754
##  -76.66226   20.38754
##  -74.38726   20.10420
##  -77.82892   20.87087
##  -78.08726   21.70420
##  -76.32059   20.16254
##  -76.63726   20.42087
##  -76.97892   20.01254
##  -78.88726   22.25420
## 
## 
## Species name:  allogus

Building an ENM

ENMTools contains functions to simplify the ENM construction process. Using enmtools.species objects and the correct modeling commands, we can build models very quickly. These commands are primarily wrappers to dismo model construction and projection functions, and at present are only available for GLM, Maxent, Domain, and Bioclim models. One of the nice bits about this setup is that it allows enmtools to automatically generate suitability maps, do model evaluation, and plot the marginal suitability of habitat for each variable separately.

GLM

GLMs usually require the user to supply a formula, an enmtools.species object, and some environmental data. If your formula is a strictly additive function of all of the environmental layers in env, though, enmtools.glm will build a formula automatically.

ahli.glm <- enmtools.glm(species = ahli, env = env, f = pres ~ layer.1 + layer.2 + layer.3 + layer.4, test.prop = 0.2)
## Adding environmental data to species ahli 
##  Processing presence points...
##  Processing background points...
allogus.glm <- enmtools.glm(species = allogus, env = env, f = pres ~ layer.1 + layer.2 + layer.3 + layer.4, test.prop = 0.2)
## Adding environmental data to species allogus 
##  Processing presence points...
##  Processing background points...
ahli.glm
## 
## 
## Formula:  presence ~ layer.1 + layer.2 + layer.3 + layer.4
## <environment: 0x7fc5b3dd05b8>
## 
## 
## Data table (top ten lines): 
## 
##       Longitude   Latitude   layer.1   layer.2   layer.3   layer.4   presence
## ---  ----------  ---------  --------  --------  --------  --------  ---------
## 1      -80.0106    21.8744      2765      1235      1174       252          1
## 4      -79.8251    21.8095      2207      1877       967       259          1
## 5      -79.8807    21.8374      2244      1828       945       249          1
## 6      -79.9550    21.8374      2250      1766       919       235          1
## 7      -80.3446    22.0136      2201      1822       978       277          1
## 8      -80.2983    21.9951      2214      1786       986       284          1
## 9      -80.1776    21.9023      2287      1722       992       266          1
## 10     -80.1591    21.9673      2984       965      1311       237          1
## 11     -80.1498    21.9858      3042       841      1371       221          1
## 12     -80.1220    21.9301      2898      1033      1231       242          1
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)])
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.4711  -0.1761  -0.1204  -0.0845   3.2446  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) 52.8490333 28.5760127   1.849   0.0644 .
## layer.1     -0.0144125  0.0070933  -2.032   0.0422 *
## layer.2     -0.0148595  0.0076379  -1.945   0.0517 .
## layer.3     -0.0006334  0.0080180  -0.079   0.9370  
## layer.4      0.0085116  0.0248357   0.343   0.7318  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 130.29  on 1011  degrees of freedom
## Residual deviance: 119.42  on 1007  degrees of freedom
## AIC: 129.42
## 
## Number of Fisher Scoring iterations: 8
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 12 
## n absences     : 1000 
## AUC            : 0.7642083 
## cor            : 0.09845587 
## max TPR+TNR at : -4.40285 
## 
## 
## Proportion of data wittheld for model fitting:  0.2
## 
## Model fit (test data):  class          : ModelEvaluation 
## n presences    : 4 
## n absences     : 1000 
## AUC            : 0.743875 
## cor            : 0.05039809 
## max TPR+TNR at : -4.517831 
## 
## 
## Suitability:  
## class       : RasterLayer 
## dimensions  : 418, 1535, 641630  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -86.90809, -74.11642, 19.80837, 23.2917  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer 
## values      : 2.691814e-06, 0.9999737  (min, max)

Notice this produces the same formula as:

ahli.glm <- enmtools.glm(species = ahli, env = env, test.prop = 0.2)
## Adding environmental data to species ahli 
##  Processing presence points...
##  Processing background points...
ahli.glm
## 
## 
## Formula:  presence ~ layer.1 + layer.2 + layer.3 + layer.4
## <environment: 0x7fc5b2c608e0>
## 
## 
## Data table (top ten lines): 
## 
##       Longitude   Latitude   layer.1   layer.2   layer.3   layer.4   presence
## ---  ----------  ---------  --------  --------  --------  --------  ---------
## 1      -80.0106    21.8744      2765      1235      1174       252          1
## 2      -79.9086    21.8095      2289      1732       957       231          1
## 3      -79.8065    21.7631      2158      1870       983       253          1
## 5      -79.8807    21.8374      2244      1828       945       249          1
## 7      -80.3446    22.0136      2201      1822       978       277          1
## 8      -80.2983    21.9951      2214      1786       986       284          1
## 9      -80.1776    21.9023      2287      1722       992       266          1
## 10     -80.1591    21.9673      2984       965      1311       237          1
## 12     -80.1220    21.9301      2898      1033      1231       242          1
## 14     -80.2148    21.9394      2329      1692      1018       269          1
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)])
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.3660  -0.1837  -0.1240  -0.0843   3.1966  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept) 59.877906  27.899368   2.146   0.0319 *
## layer.1     -0.016190   0.006936  -2.334   0.0196 *
## layer.2     -0.016814   0.007426  -2.264   0.0236 *
## layer.3     -0.002510   0.008194  -0.306   0.7593  
## layer.4      0.017902   0.023473   0.763   0.4457  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 130.29  on 1011  degrees of freedom
## Residual deviance: 120.92  on 1007  degrees of freedom
## AIC: 130.92
## 
## Number of Fisher Scoring iterations: 8
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 12 
## n absences     : 1000 
## AUC            : 0.754 
## cor            : 0.08787457 
## max TPR+TNR at : -4.56992 
## 
## 
## Proportion of data wittheld for model fitting:  0.2
## 
## Model fit (test data):  class          : ModelEvaluation 
## n presences    : 4 
## n absences     : 1000 
## AUC            : 0.74125 
## cor            : 0.05397096 
## max TPR+TNR at : -3.360277 
## 
## 
## Suitability:  
## class       : RasterLayer 
## dimensions  : 418, 1535, 641630  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -86.90809, -74.11642, 19.80837, 23.2917  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer 
## values      : 1.983417e-06, 0.9999955  (min, max)

If you want a more complicated formula, though (e.g., with interactions or polynomial effects), you’ll need to supply that manually.

ahli.glm <- enmtools.glm(species = ahli, env = env, f = pres ~ poly(layer.1, 2) + poly(layer.2, 2) + poly(layer.3, 2) + poly(layer.4, 2), test.prop = 0.2)
## Adding environmental data to species ahli 
##  Processing presence points...
##  Processing background points...
ahli.glm
## 
## 
## Formula:  presence ~ poly(layer.1, 2) + poly(layer.2, 2) + poly(layer.3, 
##     2) + poly(layer.4, 2)
## <environment: 0x7fc5ccf57238>
## 
## 
## Data table (top ten lines): 
## 
##       Longitude   Latitude   layer.1   layer.2   layer.3   layer.4   presence
## ---  ----------  ---------  --------  --------  --------  --------  ---------
## 2      -79.9086    21.8095      2289      1732       957       231          1
## 3      -79.8065    21.7631      2158      1870       983       253          1
## 4      -79.8251    21.8095      2207      1877       967       259          1
## 6      -79.9550    21.8374      2250      1766       919       235          1
## 7      -80.3446    22.0136      2201      1822       978       277          1
## 8      -80.2983    21.9951      2214      1786       986       284          1
## 10     -80.1591    21.9673      2984       965      1311       237          1
## 11     -80.1498    21.9858      3042       841      1371       221          1
## 12     -80.1220    21.9301      2898      1033      1231       242          1
## 13     -80.1776    21.9673      2914      1020      1256       237          1
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)])
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.65758  -0.16717  -0.11277  -0.07515   3.09663  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -5.0894     0.4923 -10.337   <2e-16 ***
## poly(layer.1, 2)1 -111.7489    70.4889  -1.585    0.113    
## poly(layer.1, 2)2    9.5637    28.0582   0.341    0.733    
## poly(layer.2, 2)1 -114.3667    79.0120  -1.447    0.148    
## poly(layer.2, 2)2   -5.6030    33.8869  -0.165    0.869    
## poly(layer.3, 2)1    2.0866    44.4962   0.047    0.963    
## poly(layer.3, 2)2   11.8612    17.7449   0.668    0.504    
## poly(layer.4, 2)1    3.3925    15.1844   0.223    0.823    
## poly(layer.4, 2)2   -7.8530    12.9228  -0.608    0.543    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 130.29  on 1011  degrees of freedom
## Residual deviance: 114.97  on 1003  degrees of freedom
## AIC: 132.97
## 
## Number of Fisher Scoring iterations: 8
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 12 
## n absences     : 1000 
## AUC            : 0.8186667 
## cor            : 0.118871 
## max TPR+TNR at : -4.786353 
## 
## 
## Proportion of data wittheld for model fitting:  0.2
## 
## Model fit (test data):  class          : ModelEvaluation 
## n presences    : 4 
## n absences     : 1000 
## AUC            : 0.68675 
## cor            : 0.03338417 
## max TPR+TNR at : -4.718923 
## 
## 
## Suitability:  
## class       : RasterLayer 
## dimensions  : 418, 1535, 641630  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -86.90809, -74.11642, 19.80837, 23.2917  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer 
## values      : 2.220446e-16, 0.9999816  (min, max)

To check out the marginal response functions, you only need to type

ahli.glm$response.plots
## $layer.1

## 
## $layer.2

## 
## $layer.3

## 
## $layer.4

You can also visualize your models and data in a 2D environment space using any pair of layers from your environment stack. These plots hold all non-plotted variables (layer.1 and layer.3 in this case) constant at their mean value across all presence points, then vary the plotted variables between the minimum and maximum values in env.

The suit.plot shows you suitability in environment space as a function of your two variables, with brighter colors representing variable combinations predicted to be more suitable. The points represent the occurrence points for your species in that environment space.

The colored raster of the background.plot shows you the density of background points in environment space, while the white points again represent your occurrence points in environment space.

visualize.enm(ahli.glm, env, layers = c("layer.2", "layer.4"))
## $suit.plot

## 
## $background.plot

GAM, Bioclim, Domain, and Maxent

The procedure for building Bioclim, Domain, and Maxent models is similar to the procedure for GLMs, with the exception that you do not need to pass a formula to the model function. Note that running Maxent models requires a bit of extra setup; see dismo documentation for details.

ahli.gam <- enmtools.gam(ahli, env, test.prop = 0.2)
ahli.dm <- enmtools.dm(ahli, env, test.prop = 0.2)
ahli.bc <- enmtools.bc(ahli, env, test.prop = 0.2)
ahli.mx <- enmtools.maxent(ahli, env, test.prop = 0.2)

Metrics: breadth, correlation, and overlap

ENMTools provides a number of metrics for ENMs and for similarities between ENMs. These include measures of niche breadth, based on Levins(1968). An important caveat when interpreting these metrics is that they are driven to some (variable) extent by the availability of different combinations of environmental variables. As such they are more accurately interpreted as a measurment of the smoothness of the geographic distribution of suitability scores than as an estimate of the breadth of the fundamental niche; an organism with narrow fundamental niche breadth that nonetheless encompasses a set of environmental conditions that is quite common will have a high breadth when measured using ENMs, while having a low breadth in environment space.

raster.breadth(ahli.glm)
## $B1
## [1] 0.8650396
## 
## $B2
## [1] 0.1802174

ENMTools also provides metrics for measuring similarity between ENMs. These include Schoener’s D (Schoener 1968), I (Warren et al. 2008), and the Spearman rank correlation coefficient between two rasters. While D and I are commonly used in the ENM literature, they may tend to overestimate similarity between ENMs when many grid cells are of similar values (e.g., when two species prefer different habitat but the region contains a great deal of habitat that is unsuitable for both).

raster.overlap(ahli.glm, allogus.glm)
## $D
## [1] 0.3616807
## 
## $I
## [1] 0.6267567
## 
## $rank.cor
## [1] 0.4740915

A new feature of the R version of ENMTools is that you can now use these same metrics in the n-dimensional space of all combinations of environmental variables, instead of restricting your measures of model similarity to those sets of conditions that appear in the training region. This is done by repeatedly drawing Latin hypercube samples from the space of all possible combinations of environmental variables given the min and max of each variable within the training region. ENMTools continues to draw samples until subsequent iterations differ by less than a specified tolerance value. Lower tolerance values result in more precise estimates of overlap, but can take much longer to calculate.

env.overlap(ahli.glm, allogus.glm, env, tolerance = .001)
## $env.D
## [1] 0.2753451
## 
## $env.I
## [1] 0.4720855
## 
## $env.cor
## [1] -0.00658996

Hypothesis testing

Niche identity or equivalency test

In this example, we will run a niche identity (also called equivalency) test, as in Warren et al. 2008. This test takes the presence points for a pair of species and randomly reassigns them to each species, then builds ENMs for these randomized occurrences. By doing this many times, we can estimate the probability distribution for ENM overlap between species under the null hypothesis that the two species’ occurrences in the environment are effectively a random draw from the same underlying distribution. Note that niche evolution is only one of many reasons why two species’ realized environmental distributions might cause departures from this null hypothesis. See Warren et al. 2014 for details.

To run an identity test, we need to decide what type of models we will build, how many replicates we will run, and (in the case of GLM) a model formula to use for empirical models and the Monte Carlo replicates. The resulting object contains the replicate models, p values, and plots of the results. Typically idenity tests are run with at least 99 replicates, but we are using a smaller number here for the sake of execution time.

NOTE: In order for models to be comparable, both empirical and pseudoreplicate models for the identity test are conducted with pseudoabsence points pooled for the two species being compared.

id.glm <- identity.test(species.1 = ahli, species.2 = allogus, env = env, type = "glm", nreps = 4)
id.glm
## 
## 
##  
## 
## Identity test ahli vs. allogus
## 
## Identity test p-values:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.2      0.2      0.2      0.2      0.2      0.2 
## 
## 
## Replicates:
## 
## 
##                      D           I     rank.cor       env.D       env.I      env.cor
## ----------  ----------  ----------  -----------  ----------  ----------  -----------
## empirical    0.2548982   0.5075499   -0.4020943   0.0050492   0.0299701   -0.5875723
## rep 1        0.8546444   0.9844519    0.7852184   0.7513282   0.9482161    0.8399446
## rep 2        0.7665889   0.9582296    0.3574021   0.6183091   0.8775111    0.6378208
## rep 3        0.8195256   0.9752231    0.9737111   0.7759048   0.9546163    0.9225750
## rep 4        0.9490372   0.9976710    0.9807814   0.9115978   0.9939383    0.9876159

Background or similarity test

The background or similarity test compares the overlap seen between two species’ ENMs to the overlap expected by chance if one or both species was effectively choosing habitat at random from within their broad geographic range. The purpose of this test is to correct for the availability of habitat and ask whether the observed similarity between species or populations is significantly more (or less) than expected given the available set of environments in the regions in which they occur.

NOTE: In order for models to be comparable, both empirical and pseudoreplicate models for the background test are conducted with pseudoabsence points pooled for the two species being compared.

In Warren et al. 2008, we developed this test in the context of comparing one species’ actual occurrence to the random background occurrences of the other species. This is what we call an “asymmetric” test, and in our case we did the test in both directions with the idea that we might compare the results of A vs. B background to the results of B vs. A background. This may be informative in some cases, but many people have also found this asymmetry confusing (and indeed it is often difficult to interpret). For that reason, the background test here can be conducted against a null hypothesis that is generated from “asymmetric” (species.1 vs species.2 background) or “symmetric” (species.1 background vs. species.2 background) comparisons.

Here, for instance, is a Bioclim background test using the classical asymmetric approach:

bg.bc.asym <- background.test(species.1 = ahli, species.2 = allogus, env = env, type = "bc", nreps = 4, test.type = "asymmetric" )
bg.bc.asym
## 
## 
##  
## 
## Asymmetric background test ahli vs. allogus background
## 
## background test p-values:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.2      0.2      0.2      0.2      0.2      0.2 
## 
## 
## Replicates:
## 
## 
##                      D           I    rank.cor       env.D       env.I     env.cor
## ----------  ----------  ----------  ----------  ----------  ----------  ----------
## empirical    0.1328502   0.3177390   0.0706201   0.0205455   0.1116701   0.0861670
## rep 1        0.1486335   0.3305719   0.1450342   0.0424154   0.1589587   0.0986768
## rep 2        0.1651580   0.3500385   0.1639466   0.0501244   0.1753446   0.1255899
## rep 3        0.1665392   0.3649020   0.1905156   0.0751387   0.2240976   0.1826939
## rep 4        0.1624940   0.3544069   0.1801364   0.0498466   0.1885051   0.1350149

And here is a Domain background test using the symmetric approach:

bg.dm.sym <- background.test(species.1 = ahli, species.2 = allogus, env = env, type = "dm", nreps = 4, test.type = "symmetric" )
bg.dm.sym
## 
## 
##  
## 
## Symmetric background test ahli background vs. allogus background
## 
## background test p-values:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.2      0.2      0.2      0.2      0.2      0.2 
## 
## 
## Replicates:
## 
## 
##                      D           I    rank.cor       env.D       env.I     env.cor
## ----------  ----------  ----------  ----------  ----------  ----------  ----------
## empirical    0.4929334   0.7052122   0.2916150   0.1067514   0.3106242   0.2245317
## rep 1        0.9633144   0.9983875   0.8452877   0.7789384   0.9084278   0.8025493
## rep 2        0.9683391   0.9979490   0.9703954   0.5167449   0.7288396   0.7550513
## rep 3        0.9367693   0.9887789   0.8515679   0.2534965   0.4923326   0.4736000
## rep 4        0.9212382   0.9927110   0.6998317   0.6713370   0.8399510   0.7405201

Ecospat tests

Using enmtools.species objects also provides a simplified interface to the niche equivalency and similarity tests (or identity and background tests, respectively) that were developed by Broennimann et al. (2012). These tests do not rely on ENMs, instead using kernel density smoothing to estimate density of the species in environment space. Ecospat also uses the density of the available environment to correct for availability when measuring overlaps, so that overlaps are not strictly driven by availability of combinations of environmental variables.

These tests only work with two environmental axes, so they are often done with the top two PC axes of a set of environments. In our case we’ll just pick a couple of environmental layers, though (layer.1 and layer.2). Here’s an equivalency/identity test:

esp.id <- enmtools.ecospat.id(ahli, allogus, env[[c("layer.1", "layer.2")]])
## ===========================================================================
esp.id
## 
## 
##  
## 
## Ecospat identity test ahli vs. allogus
## 
## Species    Longitude   Latitude   layer.1   layer.2
## --------  ----------  ---------  --------  --------
## ahli        -80.0106    21.8744      2765      1235
## ahli        -79.9086    21.8095      2289      1732
## ahli        -79.8065    21.7631      2158      1870
## ahli        -79.8251    21.8095      2207      1877
## ahli        -79.8807    21.8374      2244      1828
## ahli        -79.9550    21.8374      2250      1766
## 
## 
## Species    Longitude   Latitude   layer.1   layer.2
## --------  ----------  ---------  --------  --------
## ahli.bg    -80.43726   22.17087      2324      1828
## ahli.bg    -79.99559   22.02920      2660      1382
## ahli.bg    -79.89559   21.76254      2184      1849
## ahli.bg    -79.76226   21.92920      2440      1733
## ahli.bg    -80.22059   22.07087      2458      1589
## ahli.bg    -80.36226   22.01254      2240      1817
## 
## 
## Species       Longitude   Latitude   layer.1   layer.2
## -----------  ----------  ---------  --------  --------
## allogus.bg     -79.2527    22.2109      2656      1683
## allogus.bg     -78.7774    22.2241      2373      1967
## allogus.bg     -78.6189    22.2373      2317      1980
## allogus.bg     -78.1039    21.1809      2461      1669
## allogus.bg     -78.0247    21.1809      2437      1684
## allogus.bg     -77.9983    20.9301      2315      1746
## 
## 
## Species    Longitude   Latitude   layer.1   layer.2
## --------  ----------  ---------  --------  --------
## allogus    -75.28726   20.02920      1692      1774
## allogus    -75.86226   20.18754      2181      1479
## allogus    -76.66226   20.38754      2107      1720
## allogus    -74.38726   20.10420      2046      1422
## allogus    -77.82892   20.87087      2338      1703
## allogus    -78.08726   21.70420      2380      1821
## 
## 
## Species       Longitude   Latitude   layer.1   layer.2
## -----------  ----------  ---------  --------  --------
## background    -80.91226   23.27087      2495      1993
## background    -80.90392   23.27087      2504      2000
## background    -80.93726   23.26254      2500      1993
## background    -80.92892   23.26254      2500      1993
## background    -80.92059   23.26254      2509      2002
## background    -80.91226   23.26254      2504      2000
## 
## 
## ecospat.id test p-values:
##    D    I 
## 0.02 0.02

## TableGrob (3 x 2) "arrange": 6 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## 5 5 (3-3,1-1) arrange gtable[layout]
## 6 6 (3-3,2-2) arrange gtable[layout]

And here’s a symmetric background test. The difference between symmetric and asymmetric for these tests is the same as for the background tests presented above.

esp.bg.sym <- enmtools.ecospat.bg(ahli, allogus, env[[c("layer.1", "layer.3")]], test.type = "symmetric")
## ===========================================================================
esp.bg.sym
## 
## 
##  
## 
## Ecospat background test symmetric ahli vs. allogus
## 
## Species    Longitude   Latitude   layer.1   layer.3
## --------  ----------  ---------  --------  --------
## ahli        -80.0106    21.8744      2765      1174
## ahli        -79.9086    21.8095      2289       957
## ahli        -79.8065    21.7631      2158       983
## ahli        -79.8251    21.8095      2207       967
## ahli        -79.8807    21.8374      2244       945
## ahli        -79.9550    21.8374      2250       919
## 
## 
## Species    Longitude   Latitude   layer.1   layer.3
## --------  ----------  ---------  --------  --------
## ahli.bg    -80.43726   22.17087      2324       938
## ahli.bg    -79.99559   22.02920      2660      1077
## ahli.bg    -79.89559   21.76254      2184       948
## ahli.bg    -79.76226   21.92920      2440      1000
## ahli.bg    -80.22059   22.07087      2458      1028
## ahli.bg    -80.36226   22.01254      2240      1006
## 
## 
## Species       Longitude   Latitude   layer.1   layer.3
## -----------  ----------  ---------  --------  --------
## allogus.bg     -79.2527    22.2109      2656      1097
## allogus.bg     -78.7774    22.2241      2373      1067
## allogus.bg     -78.6189    22.2373      2317      1065
## allogus.bg     -78.1039    21.1809      2461       846
## allogus.bg     -78.0247    21.1809      2437       877
## allogus.bg     -77.9983    20.9301      2315       907
## 
## 
## Species    Longitude   Latitude   layer.1   layer.3
## --------  ----------  ---------  --------  --------
## allogus    -75.28726   20.02920      1692       951
## allogus    -75.86226   20.18754      2181       915
## allogus    -76.66226   20.38754      2107       907
## allogus    -74.38726   20.10420      2046       797
## allogus    -77.82892   20.87087      2338       866
## allogus    -78.08726   21.70420      2380      1002
## 
## 
## Species       Longitude   Latitude   layer.1   layer.3
## -----------  ----------  ---------  --------  --------
## background    -80.91226   23.27087      2495      1075
## background    -80.90392   23.27087      2504      1082
## background    -80.93726   23.26254      2500      1077
## background    -80.92892   23.26254      2500      1077
## background    -80.92059   23.26254      2509      1082
## background    -80.91226   23.26254      2504      1080
## 
## 
## ecospat.bg test p-values:
##    D    I 
## 0.56 0.44

## TableGrob (3 x 2) "arrange": 6 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## 5 5 (3-3,1-1) arrange gtable[layout]
## 6 6 (3-3,2-2) arrange gtable[layout]

Rangebreak tests

ENMTools also allows you to perform linear, blob, and ribbon rangebreak tests as developed in Glor and Warren 2011. The linear and blob tests are two versions of a test that permit one to ask whether the geographic regions occupied by two species are more environmentally different than expected by chance. The ribbon test, meanwhile, is designed to test whether the ranges of two species are divided by a region that is relatively unsuitable to one or both forms.

For the linear and blob tests, you call them very much like you would the identity and background tests. Here’s a linear one using GLM models:

rbl.glm <- rangebreak.linear(ahli, allogus, env, type = "glm", nreps = 4)
## 
## Building empirical models...
## Adding environmental data to species ahli 
##  Processing presence points...
##  Processing background points...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Adding environmental data to species allogus 
##  Processing presence points...
##  Processing background points...
## 
## Building replicate models...
## 
## Replicate 1 ...
## Adding environmental data to species ahli 
##  Processing presence points...
##  Processing background points...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Adding environmental data to species allogus 
##  Processing presence points...
##  Processing background points...
## 
## Replicate 2 ...
## Adding environmental data to species ahli 
##  Processing presence points...
##  Processing background points...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Adding environmental data to species allogus 
##  Processing presence points...
##  Processing background points...
## 
## Replicate 3 ...
## Adding environmental data to species ahli 
##  Processing presence points...
##  Processing background points...
## Adding environmental data to species allogus 
##  Processing presence points...
##  Processing background points...
## 
## Replicate 4 ...
## Adding environmental data to species ahli 
##  Processing presence points...
##  Processing background points...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Adding environmental data to species allogus 
##  Processing presence points...
##  Processing background points...
rbl.glm
## 
## 
##  
## 
## Linear rangebreak test ahli vs. allogus
## 
## rangebreak test p-values:
##        D        I rank.cor    env.D    env.I  env.cor 
##      1.0      1.0      0.6      0.4      0.4      0.2 
## 
## 
## Replicates:
## 
## 
##                      D           I     rank.cor       env.D       env.I      env.cor
## ----------  ----------  ----------  -----------  ----------  ----------  -----------
## empirical    0.2548982   0.5075499   -0.4020943   0.0048436   0.0298245   -0.5920514
## rep 1        0.0726801   0.2408095    0.5951234   0.5345011   0.7621649    0.5297946
## rep 2        0.2548982   0.5075499   -0.4020943   0.0046822   0.0293862   -0.5909788
## rep 3        0.1811089   0.4455980   -0.1197877   0.3190094   0.5509033   -0.0320241
## rep 4        0.2548982   0.5075499   -0.4020943   0.0050295   0.0301655   -0.5847006

And here’s a blob test using Bioclim:

rbb.bc <- rangebreak.blob(ahli, allogus, env, type = "bc", nreps = 4)
## 
## Building empirical models...
## 
## Building replicate models...
## 
## Replicate 1 ...
## 
## Replicate 2 ...
## 
## Replicate 3 ...
## 
## Replicate 4 ...
rbb.bc
## 
## 
##  
## 
## blob rangebreak test ahli vs. allogus
## 
## rangebreak test p-values:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.8      0.8      0.8      0.4      0.4      0.4 
## 
## 
## Replicates:
## 
## 
##                      D           I     rank.cor       env.D       env.I     env.cor
## ----------  ----------  ----------  -----------  ----------  ----------  ----------
## empirical    0.1328502   0.3177390    0.0706201   0.0223614   0.1171762   0.0934946
## rep 1        0.0437390   0.1507199   -0.0311702   0.1773347   0.2482751   0.2017040
## rep 2        0.1328502   0.3177390    0.0706201   0.0198827   0.1070997   0.0823804
## rep 3        0.0327940   0.1544608   -0.0172882   0.2056749   0.2907126   0.2460994
## rep 4        0.1608184   0.3470140    0.1160445   0.3570956   0.5817798   0.5617702

If you want to access the individual replicates (for instance to see how your ranges are being split up), you can find them in the list named “replicate.models” inside your rangebreak test object.

rbl.glm$replicate.models$ahli.rep.1
## 
## 
## Formula:  presence ~ layer.1 + layer.2 + layer.3 + layer.4
## <environment: 0x7fc5b38076d0>
## 
## 
## Data table (top ten lines): 
## 
##       Longitude   Latitude   layer.1   layer.2   layer.3   layer.4   presence
## ---  ----------  ---------  --------  --------  --------  --------  ---------
## 79     -74.2483    20.2963      1935      1503       923       521          1
## 64     -74.8161    20.5867      2558      1292       660       706          1
## 68     -74.6312    20.4679      2934       960       332       815          1
## 59     -74.9877    20.6528      2328      1492       779       664          1
## 61     -74.9481    20.6264      2345      1410       738       656          1
## 81     -74.1691    20.2170      1608      1722       972       445          1
## 50     -75.6611    20.9961      1828      1777      1041       497          1
## 71     -74.5916    20.4283      2944       930       364       805          1
## 78     -74.2483    20.2434      1933      1431       954       485          1
## 63     -74.8293    20.5471      2418      1234       903       604          1
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)])
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.97625  -0.00030  -0.00002  -0.00001   2.74322  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) 21.915940  10.371697   2.113  0.03460 * 
## layer.1     -0.017237   0.006085  -2.833  0.00462 **
## layer.2     -0.012974   0.005079  -2.554  0.01064 * 
## layer.3      0.005104   0.004761   1.072  0.28377   
## layer.4      0.047478   0.018312   2.593  0.00952 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 186.634  on 2015  degrees of freedom
## Residual deviance:  93.555  on 2011  degrees of freedom
## AIC: 103.55
## 
## Number of Fisher Scoring iterations: 13
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 16 
## n absences     : 2000 
## AUC            : 0.980375 
## cor            : 0.2265177 
## max TPR+TNR at : -3.739233 
## 
## 
## Proportion of data wittheld for model fitting:  0
## 
## Model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class       : RasterLayer 
## dimensions  : 418, 1535, 641630  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -86.90809, -74.11642, 19.80837, 23.2917  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer 
## values      : 2.220446e-16, 0.6496604  (min, max)

rbl.glm$replicate.models$allogus.rep.1
## 
## 
## Formula:  presence ~ layer.1 + layer.2 + layer.3 + layer.4
## <environment: 0x7fc5c82241c0>
## 
## 
## Data table (top ten lines): 
## 
##       Longitude   Latitude   layer.1   layer.2   layer.3   layer.4   presence
## ---  ----------  ---------  --------  --------  --------  --------  ---------
## 62     -74.9349    20.5339      2312      1107      1243       427          1
## 69     -74.6308    20.3725      2497      1110       792       616          1
## 54     -75.2782    20.6396      2089      1642       871       617          1
## 76     -74.2747    20.1114      1972      1243      1112       422          1
## 72     -74.4595    20.1906      2190      1196       994       502          1
## 74     -74.3803    20.1378      2060      1291       918       503          1
## 65     -74.8029    20.3359      2268      1102      1199       419          1
## 19     -78.6189    22.2373      2317      1980      1065       374          1
## 70     -74.6180    20.1378      2173      1157      1082       452          1
## 51     -75.6611    20.6792      2207      1706       875       487          1
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)])
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.3716  -0.2328  -0.1894  -0.1726   3.0655  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  2.1017704  2.4440531   0.860   0.3898   
## layer.1     -0.0011967  0.0005463  -2.190   0.0285 * 
## layer.2     -0.0015443  0.0005335  -2.894   0.0038 **
## layer.3     -0.0005267  0.0012190  -0.432   0.6657   
## layer.4     -0.0009379  0.0013137  -0.714   0.4753   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 629.56  on 3064  degrees of freedom
## Residual deviance: 618.18  on 3060  degrees of freedom
## AIC: 628.18
## 
## Number of Fisher Scoring iterations: 7
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 65 
## n absences     : 3000 
## AUC            : 0.6338564 
## cor            : 0.06187162 
## max TPR+TNR at : -3.875799 
## 
## 
## Proportion of data wittheld for model fitting:  0
## 
## Model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class       : RasterLayer 
## dimensions  : 418, 1535, 641630  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -86.90809, -74.11642, 19.80837, 23.2917  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer 
## values      : 0.006466657, 0.07945109  (min, max)

For the ribbon rangebreak test, you will need one extra thing; a third enmtools.species object representing the occurrence points (for one or both species) that fall within the ribbon of putatively unsuitable habitat. In the case of these two anoles we don’t have such a ribbon, so we’ll just simulate one based on some random points.

ribbon <- enmtools.species(species.name = "ribbon")
ribbon$presence.points <- data.frame(Longitude = runif(n = 10, min = -79, max = -78.5),
                                      Latitude = runif(n = 10, min = 21.7, max = 22.1))
plot(env[[1]])
points(ribbon$presence.points, pch = 16)

ribbon$range <- background.raster.buffer(ribbon$presence.points, 20000, mask = env)
ribbon
## 
## 
## Range raster: 
## class       : RasterLayer 
## dimensions  : 418, 1535, 641630  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -86.90809, -74.11642, 19.80837, 23.2917  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer.1 
## values      : 1, 1  (min, max)
## 
## 
## 
## Presence points (first ten only): 
## 
##  Longitude   Latitude
## ----------  ---------
##  -78.64187   21.85732
##  -78.52183   21.79001
##  -78.61741   21.71052
##  -78.68714   21.77541
##  -78.73830   22.01088
##  -78.60156   22.07664
##  -78.56221   21.87031
##  -78.88530   21.85005
##  -78.89807   21.90540
##  -78.90946   21.92327
## 
## 
## Background points not defined.
## 
## Species name:  ribbon

Now we’ll run a ribbon rangebreak test using GLM models with quadratic effects. We also need to tell it the width of the ribbons to generate for the replicates. The units for the width argument are the same units that the presence points are in; e.g., if the points are in decimal degrees you should supply the width of the barrier in decimal degrees.

rbr.glm <- rangebreak.ribbon(ahli, allogus, ribbon, env, type = "glm", f = pres ~ poly(layer.1, 2) + poly(layer.2, 2) + poly(layer.3, 2) + poly(layer.4, 2), width = 0.5, nreps = 4)
## 
## 
## No background points provided, drawing background from range raster.
## Warning in couldBeLonLat(mask): CRS is NA. Assuming it is longitude/
## latitude
## 
## Building empirical models...
## Adding environmental data to species ahli 
##  Processing presence points...
##  Processing background points...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Adding environmental data to species allogus 
##  Processing presence points...
##  Processing background points...
## Adding environmental data to species ribbon 
##  Processing presence points...
##  Processing background points...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Adding environmental data to species outside 
##  Processing presence points...
##  Processing background points...
## 
## Building replicate models...
## 
## Replicate 1 ...
## Adding environmental data to species ahli 
##  Processing presence points...
##  Processing background points...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Adding environmental data to species allogus 
##  Processing presence points...
##  Processing background points...
## Adding environmental data to species ribbon 
##  Processing presence points...
##  Processing background points...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Adding environmental data to species outside 
##  Processing presence points...
##  Processing background points...
## 
## Replicate 2 ...
## Adding environmental data to species ahli 
##  Processing presence points...
##  Processing background points...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Adding environmental data to species allogus 
##  Processing presence points...
##  Processing background points...
## Adding environmental data to species ribbon 
##  Processing presence points...
##  Processing background points...
## Adding environmental data to species outside 
##  Processing presence points...
##  Processing background points...
## 
## Replicate 3 ...
## Adding environmental data to species ahli 
##  Processing presence points...
##  Processing background points...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Adding environmental data to species allogus 
##  Processing presence points...
##  Processing background points...
## Adding environmental data to species ribbon 
##  Processing presence points...
##  Processing background points...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Adding environmental data to species outside 
##  Processing presence points...
##  Processing background points...
## 
## Replicate 4 ...
## Adding environmental data to species ahli 
##  Processing presence points...
##  Processing background points...
## Adding environmental data to species allogus 
##  Processing presence points...
##  Processing background points...
## Adding environmental data to species ribbon 
##  Processing presence points...
##  Processing background points...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Adding environmental data to species outside 
##  Processing presence points...
##  Processing background points...
rbr.glm
## 
## 
##  
## 
## ribbon rangebreak test ahli vs. allogus
## 
## rangebreak test p-values...
## 
## Species 1 vs. Species 2:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.6      0.4      1.0      0.8      0.8      1.0 
## 
## Species 1 vs. Ribbon:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.4      0.4      0.2      0.4      0.4      0.4 
## 
## Species 2 vs. Ribbon:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.2      0.2      0.2      0.2      0.2      0.8 
## 
## Outside vs. Ribbon:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.4      0.4      0.2      0.2      0.2      0.8 
## 
## 
## Replicates:
## 
## Species 1 vs. Species 2:
##                   D         I    rank.cor       env.D      env.I
## empirical 0.2742489 0.5051508  0.15455193 0.020467893 0.07339637
## rep 1     0.2723466 0.5065487  0.07376484 0.011260859 0.05189974
## rep 2     0.1141439 0.3173775 -0.01915522 0.002053111 0.01687757
## rep 3     0.2785799 0.5187816  0.01615524 0.013228309 0.06217732
## rep 4     0.4299572 0.7161948 -0.18828103 0.149242318 0.30937170
##              env.cor
## empirical -0.2556596
## rep 1     -0.4034983
## rep 2     -0.8029448
## rep 3     -0.4230139
## rep 4     -0.3242595
## 
## Species 1 vs. Ribbon:
##                    D         I  rank.cor        env.D       env.I
## empirical 0.12405130 0.2641380 0.1828724 0.0007867978 0.006713102
## rep 1     0.33688792 0.5881668 0.3995109 0.0099571507 0.041027548
## rep 2     0.22134508 0.4768959 0.5665779 0.4399150549 0.666858437
## rep 3     0.32758922 0.5735703 0.3837598 0.0257900927 0.070656140
## rep 4     0.05345405 0.1904684 0.3416105 0.0001623207 0.005240783
##              env.cor
## empirical  0.1477843
## rep 1      0.4589192
## rep 2      0.6446700
## rep 3      0.4641404
## rep 4     -0.2352111
## 
## Species 2 vs. Ribbon:
##                    D         I    rank.cor        env.D       env.I
## empirical 0.05153580 0.1803862 -0.71437687 9.308834e-05 0.003681101
## rep 1     0.45363790 0.7438772  0.27883255 5.485034e-02 0.176110109
## rep 2     0.43150795 0.7351777 -0.06515655 3.705948e-02 0.118615914
## rep 3     0.46769553 0.7511526  0.34464961 4.797471e-02 0.163117177
## rep 4     0.06043074 0.1964880  0.14811053 3.090985e-02 0.147245259
##                env.cor
## empirical -0.125132471
## rep 1     -0.280372937
## rep 2     -0.668879084
## rep 3     -0.256833051
## rep 4     -0.002399887
## 
## Outside vs. Ribbon:
##                    D         I    rank.cor        env.D       env.I
## empirical 0.06687670 0.2131445 -0.73482153 0.0001481916 0.004870576
## rep 1     0.48150293 0.7682087  0.21838188 0.0453866971 0.149897383
## rep 2     0.56320319 0.8460272  0.66302443 0.1047377506 0.245976321
## rep 3     0.48355076 0.7643736  0.29185424 0.0384509537 0.133198826
## rep 4     0.06067813 0.1981604  0.02584452 0.0350447558 0.158986867
##             env.cor
## empirical -0.140034
## rep 1     -0.472773
## rep 2     -0.420043
## rep 3     -0.451324
## rep 4     -0.069941

Note that the output table here has slope, intercept, and intercept offset.

rbr.glm$lines.df
##        slope  intercept    offset
## 1 -0.9582584 -53.211664 0.3462531
## 2 -0.2418844   2.152673 0.2572096
## 3 -1.1035993 -64.467252 0.3723180
## 4  0.9509825  94.811239 0.3449971

The intercept denotes the intercept corresponding to the CENTER of each ribbon. To get the lines denoting the edges of the ribbons (for example if you want to plot the ribbons on a map), you add and substract the offset. In other words, the top edge of the ribbon is given by y = (slope * x) + intercept + offset, while the bottom edge is given by y = (slope * x) + intercept - offset.

Literature cited

Broennimann, O., Fitzpatrick, M. C., Pearman, P. B., Petitpierre, B., Pellissier, L., Yoccoz, N. G., Thuiller, W., Fortin, M.-J., Randin, C., Zimmermann, N. E., Graham, C. H. and Guisan, A. (2012), Measuring ecological niche overlap from occurrence and spatial environmental data. Global Ecology and Biogeography, 21: 481–497. doi:10.1111/j.1466-8238.2011.00698.x

Levins, R. 1968. Evolution In Changing Environments. Monographs in Population Biology, volume 2. Princeton University Press, Princeton, New Jersey, USA.

Schoener, T. W. 1968. Anolis lizards of Bimini: resource partitioning in a complex fauna. Ecology 49:704- 726.

Warren, D.L., R.E. Glor, and M. Turelli. 2008. Environmental niche identity versus conservatism: quantitative approaches to niche evolution. Evolution 62:2868-2883. doi: 10.1111/j.1558-5646.2008.00482.x

Warren, D.L., M. Cardillo, D.F. Rosauer, and D.I. Bolnick. 2014. Mistaking geography for biology: inferring processes from species distributions. Trends in Ecology and Evolution 29 (10), 572-580. doi: 10.1016/j.tree.2014.08.003